function [postnm,postloc,postcomb,totalsamples] = ...
  TFBS_DF_Bayes_MCMC_convdiag(S,PP,offset,A0,A1,A1L,priorA,T,...
  burnin,samples,maxsamples,t,diagtype)

% [postnm,postloc,postcomb,totalsamples] = ...
%   TFBS_DF_BayesPred_MCMC_convdiag(S,PP,offset,A0,A1,A1L,priorA,T,...
%   burnin,samples,t)
%
% This function calls TFBS_DF_BayesPred_MCMC repeatedly as long as a chosen
% convergence criterion/threshold 't' for two parallel MCMC chains is met.
% 'diagtype' defines the diagnostic method: 1 = use posterior over the
% number of bindsing sites, 2 = use posterior over the binding locations.

if isempty(burnin) | isempty(samples)
  m = length(A1); % The number of motifs.
  c = 100000;
  burnin = 500000;
  samples = 500000;
end

disp(['burnin = ',num2str(burnin),', MCMC samples = ',num2str(samples)]);

[postnm1,postloc1,postcomb1,T1,a_r1] = TFBS_DF_Bayes_MCMC(S,...
  PP,offset,A0,A1,A1L,priorA,T,burnin,samples);

[postnm2,postloc2,postcomb2,T2,a_r2] = TFBS_DF_Bayes_MCMC(S,...
  PP,offset,A0,A1,A1L,priorA,T,burnin,samples);

if diagtype==1
  ct = max(abs(postnm1-postnm2));
elseif diagtype==2
  ct = max(abs(postloc1-postloc2));
else
  error('unknown ''diagtype''');
end

if ct<=t
  postnm = (postnm1+postnm2)/2;
  postloc = (postloc1+postloc2)/2;
  postcomb = (postcomb1+postcomb2)/2;
  totalsamples = samples;
else
  postnm1s = samples*postnm1;
  postloc1s = samples*postloc1;
  postcomb1s = samples*postcomb1;
  
  postnm2s = samples*postnm2;
  postloc2s = samples*postloc2;
  postcomb2s = samples*postcomb2;
  
  totalsamples = samples;
  r = 1;
  
  while 1
    disp(['Increasing MCMC iterations to ',num2str(samples+r*samples)]);
    
    [postnm1,postloc1,postcomb1,T1,a_r1] = TFBS_DF_Bayes_MCMC(S,...
      PP,offset,A0,A1,A1L,priorA,T1,0,samples);
    postnm1s = postnm1s + samples*postnm1;
    postloc1s = postloc1s + samples*postloc1;
    postcomb1s = postcomb1s + samples*postcomb1;
    
    [postnm2,postloc2,postcomb2,T2,a_r2] = TFBS_DF_Bayes_MCMC(S,...
      PP,offset,A0,A1,A1L,priorA,T2,0,samples);
    postnm2s = postnm2s + samples*postnm2;
    postloc2s = postloc2s + samples*postloc2;
    postcomb2s = postcomb2s + samples*postcomb2;
    
    totalsamples = totalsamples + samples;
    postnm1 = postnm1s/totalsamples;
    postnm2 = postnm2s/totalsamples;
    postloc1 = postloc1s/totalsamples;
    postloc2 = postloc2s/totalsamples;
    
    r = r + 1;
    
    if diagtype==1
      ct = max(abs(postnm1-postnm2));
    elseif diagtype==2
      ct = max(abs(postloc1-postloc2));
    else
      error('unknown ''diagtype''');
    end
    
    if ct<=t | totalsamples>=maxsamples
      postnm = (postnm1s+postnm2s)/(2*totalsamples);
      postloc = (postloc1s+postloc2s)/(2*totalsamples);
      postcomb = (postcomb1s+postcomb2s)/(2*totalsamples);
      postcomb1 = postcomb1s/totalsamples;
      postcomb2 = postcomb2s/totalsamples;
      break;
    end
  end
end

